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Abstract 

A newly developed spectral compressible linear 
stability code (SPECLS) using a staggered pressure 
mesh for analysis of shear flow stability is presented 
and applied to high-speed boundary layers and free 
shear flows. The formulation is the first application 
of a staggered mesh for a compressible flow anal- 
ysis by a spectral technique. For equivalent accu- 
racy of growth rates, an order of magnitude fewer 
points are needed by SPECLS than by a finite dif- 
ference formulation. Supersonic disturbances, which 
are found to have highly oscillatory structures, have 
been resolved by a spectral multi-domain discretiza- 
tion, which requires a factor of three fewer points 
than the single-domain spectral stability code. The 
study indicates, as expected, that stability of mixing 
layers is enhanced by viscosity and increasing Mach 
number. The mean flow involves a jet being injected 
into a quiescent gas. Higher temperatures of the in- 
jected gas also enhance stability characteristics of the 
free shear layer. 

Introduction 

There is an enduring interest in the linear sta- 
bility of shear flows. It can be attributed to the 
fact that currently this theory (used in conjunction 
with some semiempirical procedure such as an N- 
factor method) is virtually the only means of pre- 
dicting the location of flow transition or indeed of 
determining whether a given laminar shear flow will 
become turbulent or not. At the present time, there 
is no prospect of a unified theory of transition even in 
low-speed flows (where some of the underlying mech- 
anisms are relatively well-known), let alone in hy- 
personic flows. In view of the National Aero-Space 
Plane Project, the study of laminar-turbulent transi- 
tion in supersonic and hypersonic flows has become 
extremely important. 

As a short-term goal, it is imperative then to 
obtain linear stability results that account for fac- 
tors such as body curvature, three-dimensional mean 
states, shocks, and real gas effects. The implicit as- 
sumption is that supersonic and hypersonic transi- 
tion has its origin in linear instability and is not 
overly sensitive to details of the disturbance en- 
vironment (not that Morkovin’s bypasses (ref. 1), 
i.e., when flow transition occurs without any ini- 
tial linear instabilities, are inapplicable to the hy- 
personic regime). Clearly, as Mack (ref. 2) points 
out, efficient, accurate, and robust linear stability 
codes are needed for use in a large number of design 
calculations. 

This is the motivation behind the present work, 
which presents the first spectral compressible linear 


stability code (SPECLS) to solve the linearized 
Navier-Stokes equations. The newly developed pres- 
sure staggered mesh is shown to perform very well 
for spectral solutions to compressible flow problems. 
Spectral methods are known to yield highly accurate 
results using fewer grid points than finite difference 
formulations (refs. 3 and 4). In this study, results 
are presented to further substantiate this claim. In 
addition, a multi-domain spectral discretization in 
SPECLS (ref. 5) presented here deals economically 
with complex flows, which can include such features 
as multiple interior shear layers. The method is ver- 
ified by comparison of results for boundary layers 
from an existing finite difference compressible linear 
stability code and by comparison with analytical re- 
sults obtained for free shear flows. 

Symbols 


A 

matrix coefficients for second- 
derivative operators from equa- 
tions (14) to (18) 

B 

matrix coefficients for first- 
derivative operators from equa- 
tions (14) to (18) 

C 

matrix coefficients for zeroth- 
derivative operators from equa- 
tions (14) to (18) 

<7l,<7 2 ,<7 3 

adjustable constants for grid 
stretching function, equations (21) 
and (22) 

c* 

specific heat at constant pressure 

D 

matrix coefficient for first-derivative 
operator premultiplied by uj 

D 

differentiation with respect to y 

E 

matrix coefficient for zeroth- 
derivative operator premultiplied 
by w 

e h 

dimensional viscous stress tensor 
defined in equation (5) 

G 

quantity that determines flow 
character, defined in equation (33) 

I 

spectral interpolation operator 

L 

spectral differentiation operator 

M 

Mach number 

M c 

critical Mach number 

P 

mean flow pressure plus pressure 
perturbation 

P 

pressure perturbation 


p* 

time-dependent dimensional 

z 

spanwise coordinate 


pressure 

a 

disturbance wave number in the 

Po 

mean flow pressure 


streamwise direction 

g + 

= au + 3w 

0 

disturbance wave number in the 




spanwise direction 

<7~ 

— aw 4- (3u 

fir 

R* 

gas constant 

ratio of quiescent gas to injected gas 
free-stream temperature, T^/^l 

Re 

Reynolds number based on displace- 
ment thickness 

i 

ratio of specific heats 


dimensional time 

6* 

displacement thickness of the 

t* 


velocity profile in the x-direction 

T 

mean flow temperature plus temper- 

e 

wave propagation angle of the 


ature perturbation 


disturbance 

T e 

free-stream temperature; equivalent 
to T\ 

K* 

coefficient of thermal conductivity 


A 

second coefficient of viscosity 

To 

mean flow temperature 

P 

dimensional coefficient of viscosity 

T i 

free-stream temperature of injected 
gas 

Ve 

coefficient of kinematic viscosity 

7 2 

free-stream temperature of quies- 

P 

complex density disturbance 


cent gas 

r 

time-dependent dimensional density 

u 

streamwise mean flow plus pertur- 

Po 

mean flow density 


bation velocity 

a 

Prandtl number 

u 

streamwise complex disturbance of 

T 

complex temperature disturbance 


velocity 

ftj 

viscous stress tensor defined in 

U e 

mainstream velocity 


equation (6) 

Uo 

mean flow streamwise velocity 

<t> 

= (ui 9 + , t) T 

V 

normal mean flow plus perturbation 

G 0 

complex frequency 


velocity 

Subscripts: 


V 

Vo 

normal complex disturbance of 
velocity 

mean flow normal velocity 

G 

GL 

Gauss points 
Gauss-Lobatto points 


i 

imaginary part of 

W 

spanwise mean flow plus perturba- 
tion velocity 

0 

mean flow quantity 


spanwise complex disturbance of 

r 

real part of 

W 



velocity 

oc 

free-stream quantity 

Wo 

mean flow spanwise velocity 

Superscripts: 

X 

streamwise coordinate 

G 

Gauss points 

y 

normal coordinate 

GL 

Gauss-Lobatto points 

Vc 

computational coordinate in the 

T 

transpose 


normal direction 

* 

dimensional quantity 

l/max 

maximum extent of normal 

i 

differentiation with respect to y 


coordinate 

Acronyms: 


Vp 

normal coordinate in physical space 

COSAL 

compressible stability analysis code 
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MDSPD multi-domain spectral discretization 

SDSPD single-domain spectral discretization 

SPECLS spectral compressible linear stability 

code 

Basic Equations for Parallel Shear Flows 

Mean Flow 

The mean flows utilized in these studies are 
obtained from a spectral similarity solution for 


compressible shear flows. The basic spectral collo- 
cation technique is given in reference 6. 

For the free shear flow analysis two gases are con- 
sidered, with one of these gases taken as quiescent. 
The effect of relative temperature is studied by vary- 
ing the parameter fix = T2/T1, where T2 refers to 
the free-stream temperature of the quiescent gas and 
T\ the free-stream temperature of the injected gas. A 
typical mean flow streamwise velocity profile is given 
in figure 1. 


Compressible Stability Equations 

The basic equations governing the flow of a viscous compressible fluid are the Navier-Stokes equations. 
These equations in terms of a three-dimensional Cartesian coordinate system (xj, X2, x s) f° r a heat-conducting 
perfect gas in dimensional form are 


du[ _ t dn[ 
dt* + U ° dx* 


1 dr. 


T* 

u 


p* dx* 


d? 

dt* 


d 

dx* 


(>u*) =0 


( 1 ) 

( 2 ) 


_* * 

P c p 


dT* * dT 

+ u 


d t* 



..... .dr, 

+ T ij e 13 + Qt* dx 


b rr 


p* = p*R*T* 


where tq is the velocity in the coordinate direction x z ', 




d_r 

dx* dx* 


e‘, + (A *e* kk -p*)S, 



u 


(3) 

(4) 

(5) 

(6) 


and i,j = 1,2,3 according to the summation convention. Asterisks denote dimensional quantities and^ overbars 
denote time-dependent quantities. In these equations k* is the coefficient of thermal conductivity; R*, the gas 
constant; c*, the specific heat at constant pressure (assumed constant); p* , the first coefficient of viscosity, X , 
the second coefficient of viscosity; and the Kronecker delta. 

In this study, Xi = ( x,y,z ) and u t = (tx,v,tn). All velocities are scaled by U e , the mainstream velocity, and 
all lengths are scaled by 6*, the displacement thickness of the velocity profile in the x-direction. The Reynolds 
number and Mach number are given by 

Re = ^l (7) 


ff- t Ue (8) 

VWTe 

where u e and T e are the kinematic viscosity and mean temperature in the free stream and 7 is the ratio of 
specific heats. The Prandtl number a is assumed to be 0.70 and p* is evaluated by Sutherland’s law. 

If we assume that the base flow is a locally parallel boundary layer or free shear layer, then 


u{x,y,z,t ) = U 0 (y) + u{y)e 


t(ax+/32— ojt) 


(9) 
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( 10 ) 


v(x, y, z, t) — v(y) e t ( ax +^ z ujt ) 

w(x, y, z , I) = W 0 (y) + w(y) (n) 

p{x, y, Z, t) = P 0 (y) + p(y) ( 12 ) 

T(x, y, z, t ) = r o (y) + f(y) ( 13 ) 

where U Q , W a , P 0 , and T a represent the steady unperturbed parallel shear layer (mean flow) properties; 
quantities with tildes (~) denote complex disturbances; V a (y) is assumed zero since the flow is parallel; P a (y) 
is assumed constant across the shear layer; a and 0 are disturbance wave numbers in the x- and 2-directions, 
respectively; and u is the complex frequency. Equations (9) to (13) are substituted into equations (1) to (6), 
the mean flow terms are subtracted, and quadratic terms in the disturbances are neglected. The resulting 
system is the linearized compressible Navier-Stokes equations for the disturbance quantities given as follows 
(ref. 7): 


D 2 (au + 0w) + — < ^T' 0 D(au + 0w) + i 
i^o dl o 


■( 


2 fj, + A 
M 


l) (a 2 + 0 2 )D~v + + 0W' o )Dr 


Re 

.MoTo 


(«Uo + Wo) 


ik T ° { ° 2+s2 \ 


i Re 
Mo 


(a 2 + 0 2 )p 


zRe 

.Mo^o 


{aU 0 + j3W 0 - w ) + 


2/io “h A 
Mo . 


(a 2 + 0 2 ){au + 0w) 



dp 0 

dT 0 




aU'o) + (of/' + 0W’)T o ^ 


dTl 


(14) 


D l v + ill 


Mo 


2/i + A 


D(au + ,3w) — 


Mo 


2 Mo + A 


iRe 

(ioT 0 


(aU 0 + 0W O - UJ ) +a 2 +0 2 


+ —Dv - -Dp + i(- — — — ~-T' 0 {au + 0w) 


Mo 2/i -b A 


Mo 2 Ho 4* A / dT 0 


+ 


i dfi 0 

(2 Mo + A)/i 0 dT a 


(otU'o + (3 W' 0 )t = o 


(15) 


rpf 

Dd + i(au + 0w) ~^v + nM 2 (aU 0 + 0W O - w)p - ~(aU Q + 0W O -u)f = 0 

I o 1 0 
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(16) 



D 2 t + 2(7 - l)M 2 a aV ° 2 + ^ D(au + (3w) + ^ + 2(7 - 1 W 2 a ^-^ D(aw - fin) 


a 2 4- ft 2 


Vo dTo 


Re 


L Vo'I'o 


r' - 2t(7 - l)M 2 <j(a[/o 4- 


5 + —(7 - l)M 2 (at/ 0 + /JW„ - w)p 

Vo 

oM 2 d/i„ / , 2 


^(at/ 0 + /W 0 - w) + (a 2 + /3 2 ) - (7 - l) 2 — ^ K + ) 

^ 0 7o Vo ai 0 v 7 


1 d 2 Ho T , 2 _ 1 T ” duo 


Mo dT* 0 Mo ° <*To 


f = 0 


(17) 


D 2 (aw - Pu) + -feaW' - 0U' o )Dr + ±^T’D(aw - /?«) - ^r(«^o - 0U’ o )v 
yiodTo Vodl-o 


VqTo 


+ 


T ° ^°(aW' - PU' 0 ) + ^ — (a< - /?c/") 


Mo dT 2 


zRe 


Vo^o 


{aU 0 + pW 0 - u) + (a 2 + /3 2 ) 


(au? — /3u) — 0 


where primes and D indicate differentiation with respect to y. The equation of state is 

P , ? 

P = ~ + 

Po J- o 

where p 0 is the mean flow density and p the complex density disturbance. The boundary conditions are 


(18) 


(19) 


v, au + /3u>, oiib — /3u, f — 0 
5, au4/?u?, a xc — fiii, r — ► 0 


( y 

Vy 


y = 0, Boundary layer 


-oo, Free shear flow ) 

{y — °°) J 


(20) 


Solution Technique 

Equations (14) to (18) and (20) constitute an 
eigenvalue problem for the complex frequency u) ) once 
the disturbance wave numbers a and (3 are specified. 
Discretization of these equations in the y-direction 
forms a generalized matrix eigenvalue problem, suit- 
able for computer solution. Equations (14) to (18) 
are essentially an eighth-order system; thus the eight 
boundary conditions (eqs. (20)) are sufficient for solu- 
tion, and no boundary condition is applied to the dis- 
turbance pressure. Whatever discretization scheme is 
used must respect this arrangement, if a truly accu- 
rate solution is to be expected. Since stability anal- 
yses are inherently much more sensitive to inaccura- 
cies than a large-scale aerodynamic calculation is, an 


artificial boundary condition could potentially im- 
pact stability results. 

Staggered Discretization Scheme 

The discretization scheme used here is a spec- 
tral collocation technique, using Chebyshev polyno- 
mials as basis functions. The nodes of the variables 
5, = au 4 0w ) q~ — aw — 3ii ) and f are 

located at the Gauss-Lobatto points (the extrema of 
the last Chebyshev polynomial retained in the ex- 
pansion (ref. 4)); the energy and momentum equa- 
tions are collocated at these points. Thus discrete 
boundary conditions may be imposed for these vari- 
ables at both end points of the domain. The pressure 
nodes are located at the Gauss points (the zeros of 
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the first neglected polynomial) of a Chebyshev se- 
ries one order less than that used for the other vari- 
ables; the continuity equation is collocated at these 
points. Since no Gauss points fall on the boundary, 
we are free of any requirement of providing an artifi- 
cial numerical boundary condition for pressure, and 
we have the proper balance of number of equations 
and unknowns. 

The far-field boundary of the discretized domain 
is placed at a finite distance, typically 20 6 * to 100(5* 
from the wall or shear-layer centerline. Extensive 
sensitivity studies were performed to determine the 
effect of this finite domain truncation. 

Stretching is employed in the discretization to 
improve resolution near the wall and centerline. For 
the boundary layer, either of two stretching forms is 
used: 


or 


Vp ~ 


ymajv (G‘2 — 1 )^ l y c 

( c 2 -y‘c) Cl 


2 /ma xC^y c 

yp ~ 1 + C 3 - y c 


( 21 ) 

( 22 ) 


where y p is the coordinate in the physical space 
[0,2/max], Vc is the computational coordinate ( y c £ 
[0, 1], and Ci, C2, and C3 are adjustable constants. 
In equation (21) C2 controls the amount of stretching 
and Ci the rate. The smaller the quantity (C2 — 1), 
the stronger the stretching. In equation (22) C 3 
controls both the rate and the amount of stretching. 
The smaller C3, the stronger the stretching. In this 
work, Ci is either 4 or 6 and C2 ranges from 1.2 to 
2.0; C 3 is used between 0.01 and 0 . 03 . For the shear 
layer, equation (21) is used for stretching, yielding a 
physical space of [- y max , y max ] from y c £ [-1, 1]. 

Standard spectral collocation discretization for- 
mulas (refs. 3 and 4 ) are used to form matrix dif- 
ferentiation operators for both the Gauss-Lobatto 
(Ts f) and the Gauss ( p ) grids, with the 

selected stretching function incorporated. Mean flow 
quantities from the spectral boundary layer code of 
reference 6 are spectrally interpolated onto the new 
mesh, and derivatives of these quantities obtained 
using the differentiation operators. The generalized 
matrix eigenvalue problem which results from this 
discretization is of the form 


^GL^GL^ + ^GL^GL (<£ + Ig L p ) 

+ c GL (* + ig L p) 

= w d glL gl (4 + I§ L p) 

+ E GL (V + lgV)] ( 23 ) 


for the momentum and energy equations, and 

B C L G lg^ + C G (Ig^ + P) 

= wE G (lg L * + p) ( 24 ) 

for the continuity equation, where A, B, C, D, and 
E are matrix coefficients derived from equations ( 14 ) 
to ( 18 ), L denotes a spectral differentiation opera- 
tor, the unknown vector is <j> = [v q + q~ ?] r , 
and the vector P contains the disturbance pressure 
p. Subscripts GL and G denote location at or op- 
eration on Gauss-Lobatto and Gauss point grids, re- 
spectively; ig^ and ig^ are spectral interpolation 
matrices, from Gauss to Gauss-Lobatto points and 
vice versa. 

The unknown vectors (j> and P are collected into a 
single vector, and the matrix equations ( 23 ) and ( 24 ) 
are assembled into a large generalized matrix eigen- 
value problem for input into a standard library rou- 
tine for solution. A complex modified QZ algorithm 
(ref. 8) is used to obtain the eigenvalues of the sys- 
tem directly; this is referred to as a global search. 
The most unstable eigenvalue is then selected and 
used as an input to an inverse Rayleigh method to 
purify the eigenvalue of the effects of roundoff error 
and to obtain the solution eigenvectors. In all cases, 
the global and local (Rayleigh-iterated) eigenvalues 
agreed to better than eight decimal places. No effort 
was made to optimize the standard routines used for 
these tasks; thus, no computation times are quoted 
or compared. 

Spectral Multi-Domain Technique 

The above discretization scheme is utilized in 
the spectral multi-domain technique developed by 
Macaraeg and Streett (ref. 5 ). This technique was 
formulated to handle both advection- and diffusion- 
dominated flow problems. Extremely large differ- 
ences in discretization across an interface, through 
domain size, number of points, and stretchings, have 
been shown not to disrupt exponential-order accu- 
racy (ref. 5 ). These advantages are crucial for solving 
problems with widely disparate scales, as is the case 
for flows undergoing transition or involving chemical 
reaction (ref. 9). 

A simple one-dimensional, two-region example 
serves to illustrate the present method for interfacing 
two collocation-discretized regions. Consider the fol- 
lowing second-order, potentially nonlinear boundary- 
value problem: 


[F(U)) X - vU xx = S{U) x € [-1,1) | 

U{-l) = a U{ 1) = b J ^ 
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where JP, S, and [/ are functions of x , and subscript 
x indicates differentiation with respect to x. We wish 
to place an interface at the point x — m and have in- 
dependent collocation discretizations in the regions 
x( x ) G [— 1, m] and x^ G [m, 1]. Even though 
the point x = m is an interior point of the prob- 
lem domain, simply applying a collocation statement 
there, utilizing a combination of the discretizations 
on either side, will not work; the resulting algebraic 
system is singular, because the spectral second- 
derivative operator has two zero eigenvalues; thus 
the patching together of two spectrally discretized 
domains yields potentially four zero eigenvalues in 
the overall algebraic system. Two of these eigenval- 
ues are accounted for by imposition of boundary con- 
ditions, and one by continuity of the solution at the 
interface, but one zero eigenvalue remains in the sys- 
tem. To alleviate this difficulty, a global statement 
of flux balance is used. Rewriting equations (25) as 


[G(U)] X = S(U) (26) 


where the flux is 


G(U) = F(U) - vU x (27) 

and then integrating equation (26) from -1 to 1 
results in 


G(U)\ 


\ x—\ 


-G(U) 


+ [G] 


I X= — 1 


= f 1 S(U) 

x—m J — \ 


dx 
(28) 

If the jump in flux at the interface [G] is zero, then 
equation (28) may be written 


G(U)\ 


f m S{U ) dx = G(U) - f 1 S{U) dx 
c — — 1 J - 1 x=l Jm 

(29) 


The statement of global flux balance across the two 
regions, along with the assumption that the solu- 
tion is continuous, provides the condition necessary 
to close the equation set which results from spectral 
discretization of equations (25) in two regions. Note 
that the left side of equation (29) involves the dis- 
cretization in the region x^ 1 ) G [— l,m], while the 
right side involves the region x( 2 ) G [m, 1]. Since 
spectral collocation discretization strongly couples all 
points in their respective regions, equation (29) cou- 
ples all points in both discretizations. 

Note also that no statement is made concern- 
ing whether or not equations (25) are advection- or 
diffusion-dominated. Equations (25) are considered 
scalar equations here, although the above is extend- 
able to a system. 


A thorough treatment of the spectral multi- 
domain technique for compressible flow stability can 
be found in reference 10. 

Verification 
Boundary Layer 

For verification, calculations were performed for 
the stability analysis of compressible two-dimensional 
similarity boundary layer profiles. A spectral mean 
flow code modified for compressible flow is used for 
this purpose (ref. 6). The cases to be presented are 
for Moo = 0.00001, 4.5, and 10, assuming adiabatic 
flow. Comparison with second-mode calculations by 
Mack (ref. 2) is also presented. 

Initially a low Mach number, essentially incom- 
pressible case (Mqo — 0.00001, T e = 520° R, Re = 
2200) is analyzed. A resolution study using SPECLS 
is shown in table 1(a), which presents the eigen- 
values versus number of grid points. A similar study 
is given for COSAL, a finite difference compressible 
stability code (ref. 7), in table 1(b). Accuracy of the 
growth rate to 3 significant digits is obtained with 
45 points in the spectral code; COSAL requires ap- 
proximately 500 points for equivalent accuracy. It 
is of interest to note that converged values of growth 
rate differ in the third decimal place between the two 
codes. An independent calculation by L. M. Mack 
(Jet Propulsion Laboratory) agreed to four decimal 
places with SPECLS. Corresponding eigenfunctions 
for u are given in figure 2 for both codes. The profiles 
are virtually identical. 

A higher Mach number boundary layer profile 
perturbed by a three-dimensional disturbance is an- 
alyzed next. Conditions are M 0 0 = 4.5, T e = 520°R, 
and Re = 10 000. A resolution study for the eigen- 
value computations is given in table II for SPECLS 
and COSAL. Accuracy of growth rate to 3 signifi- 
cant digits is obtained with the spectral code using 
81 points; in table 11(b), COSAL requires approx- 
imately 800 points for equivalent accuracy. Eigen- 
functions for u are again shown for each code in 
figure 3. A multi-domain spectral discretization 
(MDSPD, ref. 5) in SPECLS with two domains 
requires one-third the number of points to obtain 
accuracy equivalent to that obtained with the 
single-domain spectral discretization (SDSPD), as il- 
lustrated in table III. The savings is significant con- 
sidering that the number of operations in SPECLS 
varies as the cube of the number of points. 

Insensitivity of these solutions to the locations of 
the far-field boundary is next illustrated with the 
MDSPD. The compressible case discussed above 
utilized an outer extent of 30<5*. The same case was 
run with this far-field boundary distance halved and 
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doubled. Results are given in table IV. A factor-of- 
four change in outer extent basically does not change 
to: the phase speed is constant to seven significant 
digits, and the growth rate is unchanged for five 
significant digits. 

First and second mode eigenfunctions for M ^ = 
10, T e — 90° R, and Re = 100000 are compared with 
compressible inviscid modes of neutral subsonic so- 
lutions (ref. 2). Plotted in figure 4 are the pressure 
eigenfunctions of the first mode from SPECLS and 
reference 2, respectively. The value of a is 1.844. 
Figure 5 displays similar plots for the second mode, 
with a = 4.9877. Nondimensionalizations and prop- 
erties differ for these comparison cases. Of relevance 
are the distinctive shapes of these p modes; that is, 
the number of zeros in p is one less than the mode 
number (ref. 2). The spectral code utilized 99 points 
for this case. A further comparison between Mack’s 
calculations and SPECLS for equivalent flow proper- 
ties and viscous flow at M 0 c =4.5 and T e = 300° R 
gave agreement to three significant figures in growth 
rate and five significant figures in phase speed. 

Shear Layer 

A further test of SPECLS is presented for a paral- 
lel shear flow. For verification, results are compared 
with those obtained by Blumen (ref. 11) for shear- 
layer instability of an inviscid compressible fluid. 

The velocity profile whose stability is investigated 
is given in dimensionless form by 

U 0 — tanhy (30a) 

W 0 = 0 (30b) 

The basic thermodynamic state is assumed con- 
stant (ref. 11). Table V displays maximum growth 
rates determined by Blumen for a range of Mach 
numbers. The corresponding growth rates from 
SPECLS (Re = 10 000) are also displayed and were 
obtained using at most 81 collocation points. To de- 
termine the maximum growth rates using SPECLS, 
a plot of growth rate versus a is displayed in figure 6 
for Mqo — 0.5 and M^c — 0.9. The peaks occur at 
the growth rates given in table V for Mqo = 0.5 and 
Moo = 0.9. Note in figure 6 that the lower Mach 
number flow is unstable to a larger range of wave 
numbers. 

It is interesting to note the effect of Reynolds 
number on the maximum growth rate. Figure 7 
displays a plot of a versus u> for = 0.9 and 

Re = 100. The maximum growth rate increases 
with increasing Reynolds number indicating that this 
shear-layer profile has an inviscid instability. 

The real and imaginary parts of the eigenfunc- 
tions u, v, and p for maximum growth rate are given 


in figure 8 for = 0.9 and Re = 10000. Each 
eigenfunction plot is normalized to one. Note that 
the real part of u displays a small structure about 
the origin. This feature disappears as the neutral 
curve is approached, as illustrated in figure 9. 

The manifestation of this structure is reflected 
in plots of magnitude and phase for u as given in 
figures 10 and 11. Note the pronounced peak in the 
magnitude which decreases radically as the structure 
of u about the inflection point disappears (fig. 10). 
As indicated in figure 11, the phase changes when the 
structure is prominent (a = 0.1) from the smooth 
phase plot when the structure has disappeared (a = 
0-4). Viscosity appears to widen the appearance of 
this structure as exhibited in figure 12, which is a 
progression of u plots for increasing a (Re = 100). 

A mean flow profile given by 

U 0 = 0.5(1 + tanhy) (31) 

is also investigated (ref. 12) and verified with the 
incompressible result of reference 9. This growth 
rate should be multiplied by 2 to compare with the 
tanh y profile growth rate since the amplitude of 
equation (31) is one-half the value of equation (30a). 
The incompressible growth rate for this profile is also 
given in table V for comparison. Eigenfunction plots 
obtained from analyzing it are given in figure 13. 

Sensitivity Analysis 

When interpreting the results of a linear stability 
analysis, one must realize the sensitivity of these re- 
sults to the accuracy of the mean flow with respect to 
the numerical scheme as well as properties assumed 
in the mean flow. To elucidate these points, we ex- 
amine the sensitivity of calculated growth rates to 
the accuracy of the mean flow and assumed thermo- 
dynamic properties. 

Consider a boundary layer flow with M 0 c = 
4.5, T e = 300°R, A = (—2/3)//, a = 0.6, and 
j3 = 1.0392. This flow defines case 1 in table VI. (A 
lower T e than in table 11(a) is chosen to emphasize 
the extreme sensitivity of results to thermodynamic 
properties even at these relatively low temperatures.) 

A less than 3-percent change in the Prandtl num- 
ber for this flow causes a 20-percent change in the 
calculated growth rate. Even more dramatic is the 
effect of changing the free-stream temperature ( T e ) 
by 42 percent, which causes a 30-percent change in 
growth rate. A lesser effect, but nonetheless impor- 
tant for code validation purposes, is the effect of bulk 
viscosity illustrated by comparing case 4 with case 1. 
The effects of these properties are of course intensi- 
fied at higher temperatures. Further points of inter- 
est are the effects of the numerical accuracy of the 
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mean flow, as well as the interpolation scheme em- 
ployed to transfer mean flow values onto the SPECLS 
grid. These latter effects, although minimal, cumu- 
latively cause a difference in the second significant 
figure of the growth rate relative to case 1. 

Also relevant is the effect of Prandtl number 
on the generalized inflection points (defined by 
{U f 0 /T o y = 0) for shear flows, since these flows 
are characterized by inflectional instabilities. For 
demonstration, we compare plots of (U' 0 /T 0 ) f versus 
y in figure 14 corresponding to two shear-layer mean 
flows with Moo ~ ^ and Pt — 0*6* The solid line 
represents a mean flow with o — 0.7 and Sutherland- 
law viscosity, representative of a typical profile in 
the present report. The dashed curve is for a — 1 
and jj, — T, representative of a profile from refer- 
ence 13. Note that the number and locations of in- 
flection points in these two cases are vastly different, 
which has important implications for studies involv- 
ing overly simplified constitutive relations. Multiple 
inflection points are found not to occur for as high 
a Mach number as 10 assuming the more realistic 
constitutive relations. 

Stability of a Compressible Shear Layer 

The spectral stability code is now used to ana- 
lyze the stability of a compressible shear layer. This 
mean flow is a spectral collocation similarity solution 
for free shear flows and is studied to understand the 
impact of transition on fuel-air mixing efficiency in 
scramjet combustors. It has been observed experi- 
mentally that the mixing efficiency is decreased four- 
fold in the range from Mach 1 to 4 (refs. 14 and 15). 
The cause of this trend is unknown. However, it is 
well-known that turbulent mixing is many orders of 
magnitude faster than laminar. Ideally one would 
like to be able to manipulate the downstream evolu- 
tion of shear layers to enhance mixing. The initial 
stages of shear-flow instabilities are driven by linear 
mechanisms. Understanding the growth and prop- 
agation of the disturbance in these early stages not 
only would allow a better understanding of the onset 
of transition but also would allow initiation of the 
transition process in a numerical model so that the 
physics might be systematically studied. This study 
investigates a range of Mach numbers, gas tempera- 
tures, and disturbance wave numbers. Unless other- 
wise stated, Re = 10000. The mean flow in all cases 
involves a jet being injected into a quiescent gas. 

Effect of Mach Number and Temperature on 
Stability 

For the cases under study, lower Mach numbers 
correspond to a more unstable mean flow than higher 


Mach numbers do. This point is illustrated in fig- 
ure 15, which displays growth rate versus a for a 
range of Mach numbers. The wave propagation an- 
gle 6 is defined as 


9 — cos 


a 

(«2 + £ 2 ) 1/2 


(32) 


The value of this angle is chosen to be 60° since 
disturbances that are propagated at an angle of 50° 
to 60° relative to the direction of flow are known 
(ref. 9) to experience the greatest amplification. As 
can be seen in figure 15, the lower the Mach number, 
the wider the band of wave numbers that can cause 
the flow to become unstable. 

The stability characteristics of the shear flow are 
also quite sensitive to the temperature difference be- 
tween the injected and the quiescent gas. Figure 16 
illustrates this sensitivity. The plot displays growth 
rate versus a for a wave propagation angle of 60° and 
Mach 3 flow. The three temperature cases given in 
the plot are cold injection (500°R) into a hot quies- 
cent gas (2500°R), fir = 5; cold injection into a cold 
gas (500°R), — 1; and hot injection (2500°R) 

into a cold quiescent gas (500° R), fix = 0.2. The 
condition relevant to the scramjet is fix = 5, since in 
this context a cold fuel is injected into a much hotter 
airstream (ref. 16). Mixing efficiency is greatly en- 
hanced by transition, so the greater instability of the 
shear layer at this temperature is a favorable scenario 
for increased fuel-air mixing. Cold injection into a 
cold quiescent gas is less unstable, and hot injection 
into a cold gas is the most stable of these cases. It 
seems that if the fuel temperature increases relative 
to the airstream, its stability is enhanced. This ob- 
servation has important implications. For example, 
cooling an aircraft by running fuel under the skin hin- 
ders mixing efficiency since heating the injected fuel 
decreases the range of wave numbers that can induce 
an instability necessary for flow transition. The im- 
pact of temperature differential on stability results 
depends heavily on the assumed temperature pro- 
files. Contrasting the above effects of temperature 
on stability with previous work (ref. 13) reveals dis- 
crepancies because of overly simplified temperature 
profiles and constitutive relations assumed in these 
earlier studies. 

Effect of Mach Number, Temperature, and 
Wave Number on Disturbance Structure 

Three-dimensional modes. It is not surprising that 
distinct differences in the shapes of temperature dis- 
turbance eigenfunctions occur as the gas tempera- 
tures are varied. To illustrate this point, figure 17 
displays temperature disturbance eigenfunctions of 
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mean flows which differ only in the temperature 
of the injected gas with respect to the quiescent 
gas. The Mach number of the injected gas is 2 and 
6 = 60°, with a = 0.0862 and 0 = 0.1493. The 
qualitative differences essentially reflect the effects 
of viscosity which varies because of the different tem- 
perature distributions for each case. 

Further, the disturbance eigenfunctions signif- 
icantly tighten in structure as a increases; thus 
greater resolution is required for the SDSPD. A pro- 
gression of eigenfunction plots for increasing a (9 = 
60°) is given in figures 18 and 19. 

Difficulties in resolution similarly occur for higher 
Mach number disturbances. The SDSPD in SPECLS 
had difficulty resolving eigenfunctions beyond 
Moo = 3.75 for 0j> = 5 and as early as M 0 0 = 3 
for p'p = 0.2. Restrictions on the allowable stretch- 
ing for single-domain spectral methods contribute to 
this difficulty. The MDSPD in SPECLS is especially 
useful in these cases. Examples of its usage for cases 
requiring very severe stretchings are given later. 

To illustrate the effect of increasing Mach num- 
ber on the three temperature cases of figure 16, plots 
of pressure disturbances are given in figure 20. The 
wave propagation angle is 60°, and the disturbances 
correspond to comparable growth rates (approxi- 
mately 10~ 2 ) in each case. Note that for fi T = 0.2 
(fig. 20), an additional lobe develops at M 0 0 = 3. Be- 
yond Moo = 3 the single-domain code has difficulties 
resolving disturbance eigenfunctions. Contrast fig- 
ure 20(a) with figure 20(b), for 0 T = 5. The p eigen- 
functions maintain a single lobe and are one-signed. 
Cases are resolved up to M^ = 3.75 (note that the 
spatial scale is for the interval [-10, 10], though the 
actual extent is [-100, 100]). Similarly, the structure 
at 0T — ^ (fig. 20(c)) is resolvable to Moo = 3.75 
although the real part of p differs from the preceding 
case. 

Two-dimensional modes . The above studies in- 
volve three-dimensional disturbances with a wave 
propagation angle of 60°. Preliminary results indi- 
cate that two-dimensional disturbances (0 = 0) ex- 
hibit similar trends; however, difficulties in resolv- 
ing eigenfunctions with the SDSPD occur at lower 
Mach numbers. Again at higher a, disturbance 
eigenfunctions tighten radically beyond the limit at 
which a single-domain spectral method can resolve, 
and cases become harder to resolve as 0 T decreases 
(cooling). This observation is demonstrated in fig- 
ure 21, which plots a sequence of p and f eigen- 
functions for increasing a. The mean flow corre- 
sponds to Moo = 1 and Re = 10000 with 0 T = 0.2. 
Note the extremely tight structures in p and t as 


q increases. Again note that while f is plotted on 
the interval [-10,10], the actual spatial extent is 
[—100, 100]. Numerical oscillations in p for a = 5 
and 6 are quite pronounced, indicating breakdown of 
resolution. 

A MDSPD is ideal for flows with radical scale 
differences (ref. 9). To illustrate this point, p and 
v eigenfunctions obtained from a multi-domain dis- 
cretization are displayed in figure 22 for Re = 10 000 
and a — 5, the case discussed previously. Adjacent 
to the MDSPD results are the SDSPD solutions for 
p and v. The multi-domain solution remains oscilla- 
tion free. The number of points in each discretization 
is roughly 100; however, the multi-domain solution 
utilizes three domains, with 41 points in the center 
domain between -0.5 and 0.5, 25 points in the left 
domain between -100 and -0.5, and 37 points in 
the right domain between 0.5 and 100. The plot of 
the pressure disturbance for the entire spatial extent 
[—100, 100] as obtained from the multi-domain solu- 
tion is given in figure 23, to illustrate the fineness of 
the structure which is resolved. 

The preceding case (a = 5) is found to be re- 
solvable by the SDSPD, but only after considerable 
trial and error stretching of the mesh for a vari- 
ety of resolutions. The important point is that the 
MDSPD is quite robust and gives accurate eigen- 
values over a wide range of stretching parameters. 
This observation is illustrated in table VII, which 
lists a range of stretching parameters and corre- 
sponding phase speeds and growth rates again for 
Mqo — 1, Re = 10000, /3j- = 0.2, and a ~ 5. Both 
the SDSPD and the MDSPD for this illustration uti- 
lized 99 points with the same outer extent. Note that 
for the SPSPD cases changing the stretching param- 
eter by approximately 25 percent causes about a 10- 
percent change in phase speed and over a 50-percent 
change in growth rate. Contrast this sensitivity with 
the MDSPD. The stretching parameter is allowed to 
vary 100 percent in the center domain (three domains 
are utilized with interfaces at dhl). The phase speed 
has changed by less than 1 percent and the growth 
rate by only 6 percent. This robustness is extremely 
important since one usually has no idea of the value 
of the phase speed or growth rate. In addition, one 
needs to determine whether the disturbance mode is 
spurious — by measuring its persistence over a wide 
range of resolutions. 

Supersonic Disturbances 

As mentioned earlier, the disturbance eigenfunc- 
tions become increasingly complex as Mach number 
is increased. These higher Mach number cases are 
unresolvable by the SDSPD. To illustrate, a series 
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of p and u eigenfunctions calculated by the SDSPD 
stability code are displayed in figures 24 and 25, 
respectively, for a disturbance wave angle of 60°. 
Note that at = 3.5 the injected gas side of 

the disturbance begins to take on an oscillatory na- 
ture; at Moo = 3.75 these oscillations are more pro- 
nounced. It is well-known that for supersonic dis- 
turbances the eigenfunction structure is oscillatory 
(ref. 12). (A supersonic disturbance occurs when the 
wave velocity of the disturbance relative to the lo- 
cal flow, in the direction of wave propagation, has 
a magnitude greater than the speed of sound.) The 
MDSPD is able to capture the structure of these su- 
personic modes with relative ease. Figure 26 displays 
two unstable supersonic modes that are associated 
with the instability of a Mach 4, /3^ = 1, free shear 
flow. The disturbance wave numbers are a = 0.30416 
and (3 = 0.2017. The MDSPD involves three do- 
mains: 105 points on the oscillatory side, 41 points 
in the inner domain, and 25 points in the outer do- 
main where the profile is smooth. Interface locations 
are ±1. Note further the level of complexity of the 
f eigenfunctions of this case. Figure 27 is a plot of 
t for both modes on a full scale from 100 to —100 
and a greatly expanded scale from 10 to —10. The 
center structure is an added complexity which also 
requires adequate resolution and further illustrates 
the necessity of a flexible discretization scheme like 
a MDSPD. 

The Mach number at which no unstable subsonic 
modes exist (i.e., all unstable modes are supersonic) 
is called the critical Mach number (M c ). Figure 28 
displays a plot of M c versus 0T for both two- and 
three-dimensional (9 = 60°) disturbances. It is 
found that Me for a three-dimensional disturbance 
at a given is equal to the two-dimensional M c 
multiplied by the reciprocal of the cosine of the 
propagation angle of the three-dimensional mode. 

An analysis from inviscid theory of compressible 
shear flows substantiates this relationship. A quan- 
tity G, defined as follows (ref. 13), determines the 
character of the flow: 


G = 


0T 


cos 


2 0 


- M d 


U >\2 


( u ° - s) 


(33) 


It can be shown (ref. 13) that the point where modes 
appear supersonic occurs at G = 0. We solve for M 
using this condition: 


M = 


cos 6 

77 - Ty' 

u ° a 


(34) 


For a two-dimensional mode 0 = 0, so that 


M = 


- (a) 2D 


The first occurrence of this relationship is the critical 
Mach number. Therefore, for a given mean flow the 
following is true: 


cos 0 

(Me) 3D _ Do ~ (^/q)3D 
(M c ) 2 d “ 1 


D 0 - (co/a) 


2D 


(35) 


It was found that at the critical Mach number for a 
given mean flow 



3D 



2D 


so that 

(A/c) 3D = (Me) 2D (cos0) 


(36) 


(37) 


Conclusions 

The first spectral collocation linear stability code 
for compressible flow (SPECLS) is presented. The 
accurate discretization (staggered pressure grid) can 
be employed in nonlinear simulations. Verification 
cases for high-speed boundary layers indicate an 
order-of-magnitude reduction in the number of points 
required to obtain equivalent accuracy in growth 
rates with a finite difference formulation. In addition, 
a multi-domain spectral discretization (MDSPD) in 
SPECLS is found to require a factor of three fewer 
points, which is significant since the operation count 
of the spectral formulation varies as the cube of the 
number of points. The highly irregular structures 
of the supersonic modes are easily resolved by a 
MDSPD, which is shown to be highly robust over a 
wide range of stretching parameters and resolutions. 

The stability analysis of a compressible shear flow 
is presented. The study indicates that for subsonic 
disturbances stability of mixing layers is enhanced 
by viscosity, increasing Mach number, and higher 
temperatures of the injected gas. The critical Mach 
number for three-dimensional disturbances is shown 
to be higher than for two-dimensional modes, the 
proportionality constant being the inverse of the 
cosine of the propagation angle. The exact value 
of this Mach number depends on the differential 
temperature between the two gases: The lower the 
value of /3 t , the lower the Mach number at which 
supersonic modes appear. 
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Table I. Temporal Eigenvalues for Low Mach Number, Incompressible Case 
[Moo = 0.00001, T e = 520°R, Re = 2200, a = 0.2, 0 = 0] 

(a) Calculated With SPECLS (SDSPD) 


Number 

of 

points 

OJy 


35 

5.972895550 x 10 -2 

4.163678340 x 10~ 3 

39 

5.982453898 x 10~ 2 

4.051750872 x 10~ 3 

41 

5.984158831 x 10 -2 

4.016799466 x 10~ 3 

45 

5.983655920 x 10~ 2 

4.022189152 x 10~ 3 

65 

5.983575125 x 10 -2 

4.023290767 x 10 -3 

95 

5.983575084 x 10~ 2 

4.023291203 x 10 -3 

200 

5.983575084 x 10" 2 

4.023291203 x 10 -3 


(b) Calculated With COSAL 


Number 

of 

points 



100 

5.983647688 x 10 -2 

3.96449684 x 10~ 3 

200 

5.981935522 x 10~ 2 

4.004225579 x 10“ 3 

300 

5.981623402 x 10 -2 

4.01151043 x 10~ 3 

500 

5.981464629 x 10~ 2 

4.015224481 x 10~ 3 

1025 

5.981396790 x 10 -2 

4.016810997 x 10~ 3 

1200“ 

5.981375400 x 10 -2 

4.017305399 x 10 -3 


“With Richardson extrapolation. 




Table II. Temporal Eigenvalues for Higher Mach Number Case 
[Moo = 4.5, T e = 520°R, Re = 10000, a - 0.6, (3 = 1.0392] 

(a) Calculated With SPECLS (SDSPD) 


Number 

of 

points 

0-V 


45 

0.49578803 

4.36603960 x 10~ 3 

51 

.49591572 

3.749838454 x 10“ 3 

65 

.496217679 

3.727396074 x 10~ 3 

81 

.496093364 

3.761833505 x 10" 3 

95 

.496106105 

3.763494422 x 10“ 3 

120 

.496104826 

3.764796852 x 10“ 3 

151 

.496104845 

3.764820857 x 10 -3 

200 

.496104846 

3.764820881 x 10“ 3 


(b) Calculated With COSAL 


Number 

of 

points 

(J r 

Vi : 

211 

0.496035465 

3.720663427 x 10~ 3 

513 

.495990048 

3.688194736 x 10 -3 

813 

.495984488 

3.684224749 x 10 -3 

1025 

.495983122 

3.683247436 x 10~ 3 

1500 

.495981891 

3.682367822 x 10“ 3 

1500 a 

.495980770 

3.681512173 x 10 -3 


a With Richardson extrapolation. 





Table III. Temporal Eigenvalues Calculated With MDSPD in SPECLS 
[Moo - 4.5, T e = 520°R, Re - 10000, a = 0.6, 0 = 1.0392] 


Number of points 



First 

domain 

Second 

domain 

Total 

u; r 


25 

17 

42 

0.496102327 

3.764737236 x 10 -3 

31 

25 

56 

.496104982 

3.764889142 x 10 -3 

37 

25 

62 

.496104984 

3.764890672 x 10 -3 


Table IV. Effect of Far-Field Boundary Location on 
Temporal Eigenvalues Calculated With MDSPD in SPECLS 

[Moo = 4.5, T € = 520°R, Re = 10 000, a = 0.6, 0 = 1.0392] 


Number of points 

Far-field 
boundary, 
6* units 


<*>i 

First 

domain 

Second 

domain 

37 

21 

15 

0.496104813 

3.764828351 x 10“ 3 

37 

31 

30 

.496104828 

3.764822175 x 10“ 3 

37 

51 

60 

.496104848 

3.764817702 x 10^ 3 


Table V. Maximum Values of 


Mach 

number 

a 

for — 

U 0 = tanh y 
(ref. 11) 

U 0 = tanh y 
(SPECLS) 

U 0 = 0.5(1 + tanh y) 
(SPECLS) 

0 

0.445 

0.190 

0.18954 

0.094688 

.1 

.433 

.187 

.18736 


.2 

.426 

.181 

.18112 


.3 

.417 

.171 

.17105 


.4 

.409 

.158 

.15760 


.5 

.397 

.141 

.14101 


.6 

.370 

.122 

.12180 


.7 

.326 

.101 

.10012 


.8 

.279 

.078 

.07760 


.9 

.208 

.055 

.05450 


1.0 

.000 

.000 

.00000 






Table VI. Sensitivity Analysis With Varying Boundary Layer Mean Flow Properties 
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Underlined entries denote property changed relative to case 1. 



Table VII. Effect of Stretching Parameter in SPECLS 
[Moo = 1.0, Re = 10 000, (3 T = 0.2, a = 5, 0 = 0] 

(a) Single-domain: 99 points on [-100, 100] 


Stretching 

parameter 

u> r 

“t 

1.1 

1.47179 

0.272236 

1.2 

1.49883 

.159461 

1.3 

1.49806 

.158226 

1.4 

1.57814 

.266635 


(b) Multi-domain: 17 points on [-100, -1], 65 points on [-1, 1], 17 points on [1, 100] 


Stretching 

parameter 

U) r 


1.6 

1.497268 

0.167511 

2.0 

1.496913 

.159498 

2.2 

1.497146 

.158781 

4.0 

1.496359 

.156986 










(a) Calculation by SPECLS (Re = 100000). 



Figure 5. Pressure fluctuation p for the second mode of neutral subsonic solution at a = 4.9877, 3 = 0, 
Moo = 10, and T e - 90°R. 
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Figure 14. Plot of the function {U^/T 0 )' versus y for two shear layer mean flows. M oc = 5 ; Pt = °- 6 - 
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Figure 17. Effect of temperature on temperature disturbance r. M Q 0 ~ 2,6 — 60°. 





















(c) /3r = 1. 


Figure 20. Concluded. 
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Single domain 


Multi-domain 



Figure 22. Disturbances v and p obtained with MDSPD and SDSPD in SPECLS for the case described 
in figure 21. a = 5. 



y 

Figure 23. Pressure disturbance p obtained with MDSPD in SPECLS displayed on entire computational 
domain for case of figure 21 for a = 5. Discretization: center domain on [-.5,5], 41 points; left domain 
on [-100, -.5], 21 points; right domain on [.5,100], 35 points. 


37 



-100 


0 

y 


100 


-100 


o 

y 


25. Normal-velocity disturbance v for increasing Mach number and with 
60°. 
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Figure 26. Pressure disturbance p for two unstable supersonic modes obtained with MDSPD in SPECLS. 
Moo = 4; 0 = 33.5°. Discretization: center domain on [-1,1], 41 points; oscillatory side, 105 points; 
smooth side, 25 points. 
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Figure 27. Temperature disturbance r for two unstable supersonic modes of figure 26 plotted on full 
domain [100,100] and expanded domain [-10,10]. 
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Pt 

Figure 28. Critical Mach number versus f) T for two-dimensional (2D) and three-dimensional (3D) 
disturbances. v 7 
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